clc;clear;close;
% 添加必要库
addpath(genpath('../utils'));

%% 参数设置
% 拍摄图
cb_folder = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\psf数据\cb\光谱图像";% 拍摄的棋盘格文件夹
noise = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\psf数据\noise1\550nm聚集位置所有中心波长图像\650nm(3.65mm).tif"; % 拍摄的噪声图像
noise2 = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\psf数据\noise1\650nm聚集位置所有中心波长图像\650nm(4.27mm).tif"; % 拍摄的噪声图像
% 其他参数
R = 17; % 模糊半径
cb_rows = 15;
cb_cols = 20;
rect = [288,197,154,105];
wave_pos = 550;% 哪个波段位置
wave_lctf = 650;% 哪个波段滤光

%% 1、读入拍摄图
% 聚焦位置数据
[imageFileNames,FileNum,WaveLength] =getImageFileNames(cb_folder);
posI = im2double(imread(imageFileNames{WaveLength == wave_pos}));
LCTFI = im2double(imread(imageFileNames{WaveLength == wave_lctf}));
% 离焦位置数据
img_noise = im2double(imread(noise));
img_noise2 = im2double(imread(noise2));

%% 2.进行图像配准
[fixedPoints, ~, ~] = detectCheckerboardPoints(posI);
[movingPoints, ~, ~] = detectCheckerboardPoints(LCTFI);
tform = fitgeotrans(movingPoints,fixedPoints,'nonreflectivesimilarity');
u_img = imwarp(img_noise2,tform,'OutputView',imref2d(size(img_noise2)));  

figure('units','normalized','outerposition',[0 0 1 1])
subplot(1,2,1)
imshow(u_img )
title('潜在清晰图像');
subplot(1,2,2)
imshow(img_noise)
title('拍摄图像');
pause(2)
close;

%% 3.进行PSF估计
k_size = 3*R;
b_img = img_noise;
u_img = u_img(rect(2):rect(2)+rect(4)-1,rect(1):rect(1)+rect(3)-1);
b_img = b_img(rect(2):rect(2)+rect(4)-1,rect(1):rect(1)+rect(3)-1);
k_init = zeros(k_size); k_init(ceil(k_size/2), ceil(k_size/2)) = 1;   
tic
psf = PSF_estimation(k_init,u_img,b_img,R);
toc

%% 4. 对PSF进行边缘零填充
psf_size = 101;
psf = padarray(psf,[(psf_size-k_size)/2,(psf_size-k_size)/2],'both');
figure('units','normalized','outerposition',[0 0 1 1])
imshow(psf,[])
exportgraphics(gcf,'PSF.png')

%% 5.进行PSF去模糊测试结果
% 后续可删的图像输出函数
imwrite(u_img,'u_img.png');
imwrite(b_img,'b_img.png');